In [2]:
# STAT 415/615 Regression (M. Baron)
# R Lab 8. Multivariate Regression
# We’ll be predicting the home sales price based on various characteristics of the home.
# For most of our analysis, we can use the same commands as in the Univariate Regression,
# but notice that the interpretation may be different.
import pandas as pd
import numpy as np
import statsmodels.api as sm
import scipy.stats as stats
import matplotlib.pyplot as plt
# Read the data
A = pd.read_csv("https://dr-baron.github.io/415-615/Data/HOME_SALES.csv")
A.head()
Out[2]:
| ID | SALES_PRICE | FINISHED_AREA | BEDROOMS | BATHROOMS | GARAGE_SIZE | YEAR_BUILT | STYLE | LOT_SIZE | AIR_CONDITIONER | POOL | QUALITY | HIGHWAY | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 360.0 | 3032 | 4 | 4 | 2 | 1972 | 1 | 22221 | YES | NO | MEDIUM | NO |
| 1 | 2 | 340.0 | 2058 | 4 | 2 | 2 | 1976 | 1 | 22912 | YES | NO | MEDIUM | NO |
| 2 | 3 | 250.0 | 1780 | 4 | 3 | 2 | 1980 | 1 | 21345 | YES | NO | MEDIUM | NO |
| 3 | 4 | 205.5 | 1638 | 4 | 2 | 2 | 1963 | 1 | 17342 | YES | NO | MEDIUM | NO |
| 4 | 5 | 275.5 | 2196 | 4 | 3 | 2 | 1968 | 3 | 21786 | YES | NO | MEDIUM | NO |
In [4]:
# Display the names of the variables
A.columns
Out[4]:
Index(['ID', 'SALES_PRICE', 'FINISHED_AREA', 'BEDROOMS', 'BATHROOMS',
'GARAGE_SIZE', 'YEAR_BUILT', 'STYLE', 'LOT_SIZE', 'AIR_CONDITIONER',
'POOL', 'QUALITY', 'HIGHWAY'],
dtype='object')
In [5]:
# Build a multivariate regression model
X = A[[
"FINISHED_AREA",
"BEDROOMS",
"BATHROOMS",
"GARAGE_SIZE",
"YEAR_BUILT"
]]
X = sm.add_constant(X)
y = A["SALES_PRICE"]
reg = sm.OLS(y, X).fit()
print(reg.summary())
OLS Regression Results
==============================================================================
Dep. Variable: SALES_PRICE R-squared: 0.736
Model: OLS Adj. R-squared: 0.733
Method: Least Squares F-statistic: 287.1
Date: Tue, 01 Sep 2026 Prob (F-statistic): 1.72e-146
Time: 18:39:51 Log-Likelihood: -2964.7
No. Observations: 522 AIC: 5941.
Df Residuals: 516 BIC: 5967.
Df Model: 5
Covariance Type: nonrobust
=================================================================================
coef std err t P>|t| [0.025 0.975]
---------------------------------------------------------------------------------
const -2962.1158 417.147 -7.101 0.000 -3781.632 -2142.600
FINISHED_AREA 0.1276 0.007 17.806 0.000 0.114 0.142
BEDROOMS -12.5494 3.894 -3.223 0.001 -20.199 -4.900
BATHROOMS 10.4191 4.945 2.107 0.036 0.704 20.135
GARAGE_SIZE 27.2366 5.930 4.593 0.000 15.587 38.886
YEAR_BUILT 1.4797 0.215 6.872 0.000 1.057 1.903
==============================================================================
Omnibus: 151.597 Durbin-Watson: 1.372
Prob(Omnibus): 0.000 Jarque-Bera (JB): 553.646
Skew: 1.300 Prob(JB): 5.99e-121
Kurtosis: 7.324 Cond. No. 4.07e+05
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.07e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
In [6]:
# What? A negative coefficient for the Bedrooms? A house with more bedrooms is cheaper?
# Answer: yes, as long as the area of the house remains constant.
# Analysis of variance table
# In R, anova(reg) gives the sequential (Type I) sums of squares.
# We can calculate the same quantities by fitting the variables sequentially.
# Model with FINISHED_AREA only
X1 = sm.add_constant(A[["FINISHED_AREA"]])
reg1 = sm.OLS(y, X1).fit()
# Model with FINISHED_AREA and BEDROOMS
X2 = sm.add_constant(A[["FINISHED_AREA", "BEDROOMS"]])
reg2 = sm.OLS(y, X2).fit()
# Model with FINISHED_AREA, BEDROOMS, and BATHROOMS
X3 = sm.add_constant(A[["FINISHED_AREA", "BEDROOMS", "BATHROOMS"]])
reg3 = sm.OLS(y, X3).fit()
# Model with FINISHED_AREA, BEDROOMS, BATHROOMS, and GARAGE_SIZE
X4 = sm.add_constant(
A[[
"FINISHED_AREA",
"BEDROOMS",
"BATHROOMS",
"GARAGE_SIZE"
]]
)
reg4 = sm.OLS(y, X4).fit()
# Full model
X5 = sm.add_constant(
A[[
"FINISHED_AREA",
"BEDROOMS",
"BATHROOMS",
"GARAGE_SIZE",
"YEAR_BUILT"
]]
)
reg5 = sm.OLS(y, X5).fit()
# Residual sum of squares for each model
RSS0 = np.sum((y - y.mean())**2)
RSS1 = np.sum(reg1.resid**2)
RSS2 = np.sum(reg2.resid**2)
RSS3 = np.sum(reg3.resid**2)
RSS4 = np.sum(reg4.resid**2)
RSS5 = np.sum(reg5.resid**2)
# Sequential sums of squares
SS_FINISHED_AREA = RSS0 - RSS1
SS_BEDROOMS = RSS1 - RSS2
SS_BATHROOMS = RSS2 - RSS3
SS_GARAGE_SIZE = RSS3 - RSS4
SS_YEAR_BUILT = RSS4 - RSS5
# Residual degrees of freedom
df_resid = reg5.df_resid
# Mean squared error
MSE = RSS5 / df_resid
# F statistics
F_FINISHED_AREA = SS_FINISHED_AREA / MSE
F_BEDROOMS = SS_BEDROOMS / MSE
F_BATHROOMS = SS_BATHROOMS / MSE
F_GARAGE_SIZE = SS_GARAGE_SIZE / MSE
F_YEAR_BUILT = SS_YEAR_BUILT / MSE
# P-values
p_FINISHED_AREA = stats.f.sf(
F_FINISHED_AREA, 1, df_resid
)
p_BEDROOMS = stats.f.sf(
F_BEDROOMS, 1, df_resid
)
p_BATHROOMS = stats.f.sf(
F_BATHROOMS, 1, df_resid
)
p_GARAGE_SIZE = stats.f.sf(
F_GARAGE_SIZE, 1, df_resid
)
p_YEAR_BUILT = stats.f.sf(
F_YEAR_BUILT, 1, df_resid
)
# Display the sequential sums of squares
anova_table = pd.DataFrame({
"Df": [1, 1, 1, 1, 1, df_resid],
"Sum Sq": [
SS_FINISHED_AREA,
SS_BEDROOMS,
SS_BATHROOMS,
SS_GARAGE_SIZE,
SS_YEAR_BUILT,
RSS5
],
"Mean Sq": [
SS_FINISHED_AREA,
SS_BEDROOMS,
SS_BATHROOMS,
SS_GARAGE_SIZE,
SS_YEAR_BUILT,
MSE
],
"F value": [
F_FINISHED_AREA,
F_BEDROOMS,
F_BATHROOMS,
F_GARAGE_SIZE,
F_YEAR_BUILT,
np.nan
],
"Pr(>F)": [
p_FINISHED_AREA,
p_BEDROOMS,
p_BATHROOMS,
p_GARAGE_SIZE,
p_YEAR_BUILT,
np.nan
]
}, index=[
"FINISHED_AREA",
"BEDROOMS",
"BATHROOMS",
"GARAGE_SIZE",
"YEAR_BUILT",
"Residuals"
])
anova_table
Out[6]:
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| FINISHED_AREA | 1.0 | 6.655486e+06 | 6.655486e+06 | 1310.621461 | 9.450540e-144 |
| BEDROOMS | 1.0 | 2.761256e+04 | 2.761256e+04 | 5.437562 | 2.009207e-02 |
| BATHROOMS | 1.0 | 1.427102e+05 | 1.427102e+05 | 28.102991 | 1.708199e-07 |
| GARAGE_SIZE | 1.0 | 2.249873e+05 | 2.249873e+05 | 44.305289 | 7.197206e-11 |
| YEAR_BUILT | 1.0 | 2.398082e+05 | 2.398082e+05 | 47.223868 | 1.831694e-11 |
| Residuals | 516.0 | 2.620307e+06 | 5.078115e+03 | NaN | NaN |
In [7]:
# FINISHED_AREA alone explains 6655486.
# BEDROOMS explains an additional amount of 27613. Etc.
In [8]:
# Confidence intervals for the slopes.
CI = reg.conf_int(alpha=0.10)
CI.columns = ["5 %", "95 %"]
CI
Out[8]:
| 5 % | 95 % | |
|---|---|---|
| const | -3649.496034 | -2274.735649 |
| FINISHED_AREA | 0.115787 | 0.139404 |
| BEDROOMS | -18.965525 | -6.133335 |
| BATHROOMS | 2.270205 | 18.568062 |
| GARAGE_SIZE | 17.465418 | 37.007864 |
| YEAR_BUILT | 1.124881 | 1.834506 |
In [9]:
# Confidence intervals for the slopes with Bonferroni adjustment
# (just 5 slopes; suppose we are not interested in the interval for the intercept).
# Bonferroni adjustment: alpha = 0.10 / 5
CI_Bonf = reg.conf_int(alpha=0.10 / 5)
CI_Bonf.columns = ["1 %", "99 %"]
CI_Bonf
Out[9]:
| 1 % | 99 % | |
|---|---|---|
| const | -3935.569039 | -1988.662644 |
| FINISHED_AREA | 0.110873 | 0.144318 |
| BEDROOMS | -21.635767 | -3.463093 |
| BATHROOMS | -1.121206 | 21.959473 |
| GARAGE_SIZE | 13.398843 | 41.074440 |
| YEAR_BUILT | 0.977216 | 1.982171 |
In [10]:
# Testing several slopes in one hypothesis.
# H0: β4 = 0 and β5 = 0 versus H1: either β4 ≠0 or β5 ≠0
# Consider a reduced model without these variables. Compare two models via a partial F-test.
# Reduced model without GARAGE_SIZE and YEAR_BUILT
X_reduced = sm.add_constant(
A[[
"FINISHED_AREA",
"BEDROOMS",
"BATHROOMS"
]]
)
reg_reduced = sm.OLS(y, X_reduced).fit()
# Residual sum of squares for the reduced and full models
RSS_reduced = np.sum(reg_reduced.resid**2)
RSS_full = np.sum(reg.resid**2)
# Number of restrictions
q = 2
# Residual degrees of freedom for the full model
df_full = reg.df_resid
# Partial F statistic
F_partial = (
((RSS_reduced - RSS_full) / q)
/ (RSS_full / df_full)
)
F_partial
Out[10]:
45.76457820587542
In [11]:
# P-value for the partial F-test
p_partial = stats.f.sf(F_partial, q, df_full)
p_partial
Out[11]:
5.050154701448704e-19
In [12]:
# Residual plots
fig = plt.figure(figsize=(10, 8))
sm.graphics.plot_regress_exog(
reg,
"FINISHED_AREA",
fig=fig
)
plt.show()
